C...................... VPERP.FOR ................................
C... This subroutine provides the ExB drift of flux tubes in the 
C... equatorial plane  either from the teardrop model or by reading
C... by reading Vperp from a file. The data file is read in a 2-D 
C... array MEDATA(N,M) and then interpolated to model time
C... JUNIT=23 is the unit number of the data file.
C... Dr. P. Richards, Y long March 1998
C... Modified in November 2008 to read Vperp instead of L from the file.
C... Conversion of vperp to L shell is now done in theis program
      SUBROUTINE
     >   GET_VPERP(NWTIM,  !.. number of times in file (<0 on 1st call)
     >                DT,  !.. time step in seconds
     >             UTHRS,  !.. universal time for drift determination
     >               SAT,  !.. input local time for drift
     >            RUNHRS,  !.. length of the run (hours).
     >               PCO,  !.. L value to return
     >              BLON,  !.. Magnetic longitude NO LONGER returned
     >            L_STAG,  !.. stagnation point for teardrop ExB model
     >           L_SHORT,  !.. L on the short axis for teardrop model
     >            L_PHASE, !.. Time of stagnation for teardrop model
     >            VPERPN,  !.. perp velocity to return (m/s)
     >             NCOLS,  !.. # of columns read from data file.
     >            UTORLT,  !.. tells if time is UT or LT. 
     >           PCO_MIN,  !.. minimum L value with teardrop model
     >           PCO_MAX,  !.. maximum L value with teardrop model
     >           PCO_LIM)  !.. maximum L value with teardrop model

      IMPLICIT NONE
      INTEGER NWTIM, NCOLS,UTORLT,N,M,JUNIT,IPT,I,J
      INTEGER NVTIM          !.. number of interpolated ExB times
      PARAMETER (N=59999, M=4, JUNIT=30)
      DOUBLE PRECISION VRESULT(M),RUNHRS,VPERPN,PCO,PCO_SAV,DT,PCO_MIN
      DOUBLE PRECISION PCO_MAX,PCO_LIM
      REAL MEDATA(N,M),SAT,UTHRS,BLON,L_STAG,L_SHORT,SAT_DT,DEL_TIM
     >  ,L_PHASE
      REAL UTSTOP           !.. specify UT to end run in hours
      REAL LSTOP            !.. specify L-value at end of run at UTSTOP
      REAL VTIME(N),VEXB(N) !.. UT and ExB velocity from file
      REAL LVSCAL           !.. For scaling ExB drift with L-shell
      SAVE MEDATA

      DATA PCO_SAV/0.0/
      VPERPN=0.0

      !... use teardrop model Khazanov et al. JGR 1994, page 5923
      IF(L_STAG.GT.1.1) THEN
          IF(L_STAG.LT.2*L_SHORT) THEN
             WRITE(6,*) ' ** In ExB model L_SHORT < 0.5*L_STAG. Check'
             WRITE(6,*) ' ** GLAT/L-shell and L_STAG in Fyyyyddd.* file'
             CALL RUN_ERROR  !.. print ERROR warning in output files
             STOP
          ENDIF
          !.. get maximum L at dusk for setting up proper altitude grid
          CALL TEARDROP(L_STAG,L_SHORT,18.0,PCO_MAX,VPERPN)
          !.. get maximum L at dusk for setting up proper altitude grid
          CALL TEARDROP(L_STAG,L_SHORT,6.0,PCO_MIN,VPERPN)
          !.. Two calls to ExB model to evaluate Vperp = DL/DT
          SAT_DT=SAT-DT/3600.0
          CALL TEARDROP(L_STAG,L_SHORT,SAT_DT+L_PHASE,PCO_SAV,VPERPN)
          CALL TEARDROP(L_STAG,L_SHORT,SAT+L_PHASE,PCO,VPERPN)
          IF(PCO.LT.PCO_LIM) PCO=PCO_LIM   !.. Stop L-value from getting too small
          IF(PCO.LT.PCO_LIM) VPERPN=0.0    !.. If L-value does not change
          !... VPERPN is the instantaneous velocity at L. Better to use
          !... an average over the time step DT
          IF(DT.GT.1) VPERPN=(PCO-PCO_SAV)*6.37E8/DT
          NCOLS=-1
          RETURN
      ENDIF

      !..If this is the first call to GET_VPERP, call READ_VPERP to read
      !.. the MEDATA and NWTIM from the data file
      IF (NWTIM .LT. 0) THEN
        CALL READ_VPERP(N,M,JUNIT,SAT,UTHRS,RUNHRS,NVTIM,NCOLS,
     >                 UTORLT,VTIME,VEXB,LVSCAL)
        IF(NCOLS.LT.2) THEN
           WRITE(6,*) '   ** ERROR-ERROR: ExB file is bad or missing'
           CALL RUN_ERROR  !.. print ERROR warning in output files
           STOP
        ENDIF
        UTSTOP=UTHRS+RUNHRS  !.. Specify the time to stop 
        LSTOP=PCO            !.. The input L-shell is used as the final L-value
        !.. Use the vperp read in to set up the L values array for the FLIP run
        PCO_MIN=9999
        CALL Trace(N,M,NVTIM,UTSTOP,LSTOP,VTIME,VEXB,NWTIM,PCO_MAX,
     >     PCO_MIN,PCO_LIM,MEDATA,LVSCAL)
      ENDIF
      
      !.. If column number is < 2, no data file
      IF(NCOLS.NE.2) NWTIM=0
      IF(NCOLS.NE.2 .OR. NWTIM.EQ.0) RETURN

      !.... If UTORLT =1, UT is used in the hemisphere else if UTORLT=0
      !.... calculate the winds at the specified local time SAT.
         IF (UTORLT .EQ. 1) THEN
            CALL INT_VPERP(N,M,NWTIM,NCOLS,UTHRS,MEDATA,VRESULT,IPT)
         ELSE IF (UTORLT .EQ. 0) THEN
            CALL INT_VPERP(N,M,NWTIM,NCOLS,SAT,MEDATA,VRESULT,IPT)
         ELSE
            WRITE(6,*) 'Problem with UT/LT switch value in data file'
            CALL RUN_ERROR  !.. print ERROR warning in output files
            STOP
         ENDIF

      PCO=VRESULT(2)   !.. L-shell to return
      IF(PCO.LT.PCO_LIM) PCO=PCO_LIM   !.. Stop L-value from getting too small

      !.. Vperp is calculated from dL/dt
      DEL_TIM=3600.0*(MEDATA(IPT+1,1)-MEDATA(IPT,1))
      !... drift from values from file
      IF(DEL_TIM.GT.0)
     >  VPERPN=(MEDATA(IPT+1,2)-MEDATA(IPT,2))*6.37E8/DEL_TIM
      !... calculate drift from actual change in L after 1st time step
      IF(PCO_SAV.GT.0) VPERPN=6.37E8*(PCO-PCO_SAV)/DT
      IF(PCO.LT.PCO_LIM) VPERPN=0.0    !.. If L-value does not change

      PCO_SAV=PCO  !... save L value for vperp calculation

      RETURN
      END
C
C::::::::::::::::::::::: SUBROUTINE READ_VPERP :::::::::::::::::::::::::::::
C.... This subroutine reads the UT/LT and the meridional component of the E x B 
C.... Drift from a file 
C.... Three sentinel numbers are: 
C....  1. UTORLT which is used to decide UT/LT;
C....  2. NCOLS tells the # of columns in the data file. 
C.... Modified in December 2008 to read vperp rather than L-shell from data file
      SUBROUTINE 
     >   READ_VPERP(N,   !.. # of columns dimension
     >              M,   !.. # of rows dimension
     >          JUNIT,   !.. input data file unit
     >            SAT,   !.. LT in hours
     >          UTHRS,   !.. UT in hours
     >         RUNHRS,   !.. Length of run in hours 
     >          NWTIM,   !.. # of times
     >          NCOLS,   !.. # of colums in file
     >         UTORLT,   !.. UT or LT indicator
     >          VTIME,   !.. Time from data file
     >           VEXB,   !.. Velocity from data file
     >         LVSCAL)   !.. ExB drop off with L-shell

      IMPLICIT NONE
      INTEGER NCOLS,UTORLT,NWTIM,JUNIT,IOPEN
      INTEGER N,M,I,K
      DOUBLE PRECISION RUNHRS   !.. Length of run in hours
      REAL UTHRS                !.. UT in hours
      REAL HEADR                !.. for reading sentinel
      REAL VTIME(N),VEXB(N)     !.. UT and ExB velocity from file
      REAL VFACT                !.. UT and ExB velocity from file
      REAL SAT                  !.. Local time (hours)
      REAL TTIME(N)             !.. Temporary UT time (hours)
      REAL TVEXB(N)             !.. Temporary ExB velocity from file
      REAL DAYADD               !.. for adding days to TTIME
      REAL LVSCAL               !.. For scaling ExB drift with L-shell
      CHARACTER *80 CHEADR      !.. string for header

      !... Test to see if the E x B data file exists
      CALL TFILE(JUNIT,IOPEN)
      IF(IOPEN.EQ.0) THEN
         NCOLS=0
         RETURN
      ENDIF

      !.. Read first header and junk it. Must be at least 1 header 
      READ(JUNIT,'(80A)') CHEADR
       
      !.. Read remaining input file headers until sentinel value.
      !.. keep looping on error until valid sentinel is read 
 10   READ(JUNIT,'(F8.0)', ERR=10,END=20) HEADR 

 20   CONTINUE

      UTORLT=NINT(HEADR)   !.. now convert input to integer

      !.. Read other control parameters
      READ(JUNIT,*) NCOLS           !.. number of columns in the data file
      READ(JUNIT,*) VFACT           !.. velocity scaling factor
      !.. ExB variation with L-shell: V=Vperp*EXP(LSCAL*L-shell)
      READ(JUNIT,*) LVSCAL

      IF(NCOLS.NE.2) CLOSE(JUNIT)
      IF(NCOLS.NE.2) RETURN

      READ(JUNIT,'(80A)') CHEADR    !.. Read column header and junk it.

      !.. Read in the TIMES and Vperp DATA and test for bad EXB data
      DO  25 I = 1, N 
        READ(JUNIT, *, END = 26, ERR = 26) VTIME(I),VEXB(I)
        VEXB(I)=VEXB(I)*VFACT
        !.. Test Vperp is within the correct range.
        IF ((VEXB(I).GT.100000.0).OR.(VEXB(I).LT.-100000.0)) THEN
           WRITE(6,100) VEXB(I),I
 100       FORMAT(/1X,'ERROR! v=ExB is ',1P,E10.2,' OUT OF RANGE!' 
     >       ,/1X,'** Check the data file at row ',I3)
           CALL RUN_ERROR  !.. print ERROR warning in output files
           STOP
        ENDIF
 25   CONTINUE
 26   CONTINUE


      IF(I.GE.N) THEN
         CALL RUN_ERROR  !.. print ERROR warning in output files
         WRITE(6,100) 
     >      '  *** too many values in ExB drift file: unit # =',JUNIT
         STOP
      ENDIF

      NWTIM = I-1

      !... Test for the times to be in ascending order
      DO 30 I = 2,NWTIM
        IF (VTIME(I) .LT. VTIME(I-1)) THEN
          WRITE(6,781) JUNIT,VTIME(I),VEXB(I)
 781      FORMAT(//3X,'ERROR! TIME NOT IN ASCENDING ORDER in ExB'
     >      ,/3X,'data file unit=',I3,' near the following line'
     >      ,/3X,F7.2,1P,9E10.2/)
          CALL RUN_ERROR    !.. print ERROR warning in output files
          STOP
        ENDIF 
30    CONTINUE

      !.. Test for UTHRS larger than the first TIME
      IF ((VTIME(1).GT.UTHRS).AND.(IABS(UTORLT).EQ.1)) THEN
         WRITE (6,391) JUNIT,UTHRS,VTIME(1) 
 391     FORMAT(2X,'ERROR! Need v=ExB at earlier times in data'
     >     ,1X,'file unit=',I3,/5X
     >     ,'Run starts at ',F8.2,' hrs but first data time is ',F8.2)
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP
      ENDIF

      !.... Test for the total hours of run (UTHRS + RUNHRS) less
      !.... than the last TIME
       IF ((VTIME(NWTIM).LE.(UTHRS+RUNHRS)) .AND.
     >    (IABS(UTORLT).EQ.1)) THEN
         RUNHRS=VTIME(NWTIM)-UTHRS-0.02
         WRITE (6,392) RUNHRS
 392     FORMAT(/2X,'***WARNING***  run time exceeds the last time in'
     >    ,1X,' ExB file. TSTOP has been reset to ',F10.2,' hours'/)
         IF(RUNHRS.LE.0.0) THEN
            WRITE(6,*) ' NOT ENOUGH VALID DATA'
            CALL RUN_ERROR  !.. print ERROR warning in output files
            STOP
         ENDIF
       ENDIF

      !.... Convert local time to universal time. Also make sure local times are OK
      IF (IABS(UTORLT) .EQ. 0) THEN
         !.. Make sure local times are OK
         IF (ABS(VTIME(1)) .GT. 0.001)  WRITE(6,393)
 393     FORMAT(2X,'ERROR! First local time must be zero')
         IF (ABS(VTIME(NWTIM)-24.0) .GT. 0.001) WRITE(6,394)
 394     FORMAT(2X,'ERROR! last local time must be 24.0')
         IF (ABS(VTIME(1)).GT.0.001 .OR.
     >      ABS(VTIME(NWTIM)-24.0).GT.0.001) STOP
         VTIME(1) = 0.0       !.. Make sure first UT to allow any start time
         VTIME(NWTIM) = 24.0

         !.. Now Convert LT to UT. Use temporary arrays to store 
         DAYADD=-24.0    !.. DAYADD used to make UT continuous
         DO I=1,999999
           K=MOD(I-1,24) +1
           IF(K.EQ.1) DAYADD=DAYADD+24
           TTIME(I)=VTIME(K)+DAYADD+UTHRS-SAT
           TVEXB(I)=VEXB(K)
           IF(TTIME(I).GT.UTHRS+RUNHRS+48.0) GO TO 40
         ENDDO
 40      CONTINUE
         NWTIM=I-1     !.. Readjust number of times

         !.. Transfer data back to actual arrays and make UT begin at 0.0 
         !.. by subtracting 24 hours
         K=0
         DO I=1,NWTIM
           IF(TTIME(I)-24.0.LT.0) K=0  !.. Wait until time is positive
           K=K+1
           VTIME(K)=TTIME(I)-24.0
           VEXB(K)=TVEXB(I)
         ENDDO
         UTORLT=1.0    !.. indicate LT has been converted to UT
         VTIME(1)=0.0  !.. Make sure first UT is 0.0 to allow any start time
      ENDIF
      CLOSE(JUNIT)
      RETURN
      END
C:::::::::::::::::::: TEARDROP ::::::::::::::::::::::::::::
C.....This subroutine calculates L-values based on the eq.(4) of the
C.....paper by Khazanov et al JGR, page 5924, 1994. It also computes
C.....the drift velocity
C.....Written for Dr. Richards by Jinhua Liao, August 15, 1998
      SUBROUTINE 
     >   TEARDROP(Lst,  !.. equipotential location on dusk axis
     >             LN,  !.. equipotential location on noon meridian
     >             LT,  !.. Local Time at equatorial plane
     >              L,  !.. Calculated L value
     >              V)  !.. drift velocity (m/s)

      IMPLICIT NONE
      DOUBLE PRECISION L,V
      REAL  Lst, LN, LT,RE,PI,DELTA,LT24
      REAL  LAMDA        !.. longitudinal angle measured from midnight 

      LT24=AMOD(LT,24.0) !.. make sure LT is less than 24
      Re=6370000.        !..Earth radii in meter
      PI=4*ATAN(1.0)
      Lamda=2.*LT24*PI/24.
      L=LN               !.. default for noon and midnight cases

      !... If not noon nor midnight
      IF(AMOD(LT24,12.0).NE.0.0) THEN
        DELTA=Lst/(2.*LN)
        !..Test if the number under square root is positive
        IF((DELTA**2+SIN(LAMDA)).LT.0.) THEN
          WRITE(6,*) ' ExB problem: initial L is too big'
          STOP
          CALL RUN_ERROR  !.. print ERROR warning in output files
        ELSE
          !..Calculate L value based on eq.(4)
          L=Lst*(SQRT(DELTA*DELTA+SIN(LAMDA))-DELTA)/sin(LAMDA)
        ENDIF
      ENDIF
    
      !.. Calculate drift velocity making sure demoninator not zero  
      IF(2*L*SIN(LAMDA)+Lst**2/LN.NE.0.0)THEN
         V=-Re*L**2*COS(LAMDA)*PI/(2*L*SIN(LAMDA)+Lst**2/LN)/4.32E4
      Else 
         V=0.
      ENDIF

      V=V*100.0    !.. convert to cm/sec

      RETURN
      END	      
C::::::::::::::::::::::::::: INT_VPERP :::::::::::::::::::::::::::::::::::
C.... This subroutine calculates the values of VRESULT (including wind,
C.... density, temparture etc..) at required time TIME by using the 
C.... subroutine FTW. It calls BISPLT to locate the data used for 
C.... interpolation. MEDATA(PDT,NCOLS) and MEDATA(PDT+1, NCOLS) will
C.... be used for interpolation. 
C.... Dr. Phlip Richards, Kevin Zou  at July 8,1993
      SUBROUTINE INT_VPERP(N,M,NWTIM,NCOLS,TIME,MEDATA,VRESULT,PDT)

      INTEGER PDT,NWTIM,BEGIN
      DOUBLE PRECISION VRESULT(M)
      REAL MEDATA(N,M),TIME

      !..Use subroutine BISPLW to locate the data used for interpolation.
      BEGIN=1
      CALL BISPLW(N,M,BEGIN,NWTIM,TIME,MEDATA,PDT)
      !.... Use subroutine FTW to interpolate at TIME for VRESULT
      CALL FTW(N,M,NWTIM,NCOLS,TIME,PDT,MEDATA,VRESULT)
      RETURN
      END

C:::::::::::::::::::::::: LNVEL :::::::::::::::::::::::
C.....This subroutine calculates Ln based on L value and LT(local time).
C.....Given L value at a particular time, Ln can be calculated from eq.(4)
C.....of the paper by G.V.Khazanov et al(1994). Further, the drift velocity 
C.....was calculated.
C.....Written by Dr. Richards, Jinhua Liao, August 15, 1998
      SUBROUTINE 
     >     LNVEL(Lst,     !.. The location along the dusk meridian
     >           LT,      !.. Local Time
     >           L,       !.. Geocentric distance expressed in Earth radii
     >           Ln,      !.. The intersection of the equipotential 
     >                    !.. surface with the noon-midnight meridian
     >           V)       !.. drift velocity (m/s)

      REAL  Lst, LT, Ln
      REAL  LAMDA        !.. longitudinal angle measured from midnight 
      DOUBLE PRECISION L,V

      Re=6370.*1000.     !..Earth radii in meter
      TF=3600.           !..For transfering time unit from hour to second 
      PI=4*ATAN(1.0)
      Lamda=2.*LT*PI/24.
      !... Calculates Ln based on eq.(4)
      Ln=-Lst**2*L/(L**2*sin(LAMDA)-Lst**2)
      !.. Start to calculate drift velocity
      !.. Make sure that the demoninator in the velocity equation is not zero  
      IF(2*L*SIN(LAMDA)+Lst**2/Ln.NE.0.0)THEN
         V=-Re*L**2*COS(LAMDA)*PI/(2*L*SIN(LAMDA)+Lst**2/Ln)/(12.*TF)
      Else 
         V=0.
      ENDIF
      RETURN
      END	
C:::::::::::::::::::: GET_PHASE ::::::::::::::::
C... This function determines the phase of the tear drop model
C... Based on the Ap/Kp
      REAL FUNCTION TD_PHASE(UTHRS,AP2,AP3)
      IMPLICIT NONE
      REAL UTHRS,DELTIM,AP2,AP3,KP2,KP3,KPAVE
      !... determine time elapsed since last Kp/Ap change which
      !... happens every 3 hours
      DELTIM=AMOD(UTHRS,3.0)
      KP2=0
      KP3=0
      IF(AP2.GT.2) KP2=INT(ALOG(AP2)*1.7-0.84)
      IF(AP3.GT.2) KP3=INT(ALOG(AP3)*1.7-0.84)
      KPAVE=((KP2*DELTIM+(3-DELTIM)*KP3))/3.0
      TD_PHASE=42.0-(11.3+47.0/(KPAVE+3.9))
      WRITE(6,'(9F9.2)') DELTIM,KP2,KP3,KPAVE,TD_PHASE
      RETURN
      END
C..::::::::::::::::::::::::::: Trace::::::::::::::::::::::::::::::::::::::::::::
C.. FUNCTION    : Trace
C.. DESCRIPTION : This function is used to trace back the ExB drift and 
C.. find its start locations for given UT and L values. 
C.. AUTHOR      : Jinhua Liao in January, 1999, Modified for FLIP by 
C.. P. Richards Dcember 2008
C..***********************************************************
      SUBROUTINE 
     >Trace(NDIM,   !.. # of ExB times
     >      MDIM,   !.. # of columns dimension
     >     NWTIM,   !.. # of rows dimension
     >    UT_end,   !.. specified Universal Time
     >     L_end,   !.. specified L_shell
     >        UT,   !.. array of Universal Times
     >      Vion,   !.. array of ion velocity
     >     NDint,   !.. Number of data points in array (output)
     >   PCO_MAX,   !.. maximum L value
     >   PCO_MIN,   !.. minimum L value
     >   PCO_LIM,   !.. L value limit
     >    MEDATA,   !.. output for Vperp times, and velocities
     >    LVSCAL)   !.. ExB drop off with L-shell
     
      IMPLICIT NONE 
      INTEGER i,k,j,kmin       !.. loop control variables
      INTEGER NDint            !.. # of interpolated values
      INTEGER NDIM,MDIM        !.. Dimensions of MEDATA
      INTEGER NWTIM            !.. # of ExB times
      REAL time_step           !.. time step for the L values
      REAL Vint                !.. integrated drift velocity
      REAL UT_end              !.. specified Universal Time
      REAL L_end               !.. specified L_shell
      REAL LSCAL,LVSCAL        !.. Scaling factors to reduce ExB drift with L
      REAL UT(NDIM)            !.. array of Universal Times
      REAL Vion(NDIM)          !.. array of ion velocity
      REAL Veq(NDIM)           !.. array of velocity at equtorial plane
      REAL L(NDIM)             !.. array to store the start location
      REAL time(NDIM)          !.. time for specifying L
      REAL FTAE                !.. Interpolation function
      REAL Vel_Convert         !.. Function to scale velocity
      REAL Vsave(NDIM)         !.. Save the input velocity for testing
      REAL MEDATA(NDIM,MDIM)   !.. Data array
      DOUBLE PRECISION PCO_MAX !.. L value maximum
      DOUBLE PRECISION PCO_MIN !.. L value minimum
      DOUBLE PRECISION PCO_LIM !.. L value lower limit

      DATA vint/0/,time_step/0.02/


      k=NDIM              !.. intialize time loop variable to maximum                  
      L(k) = L_end        !.. assign the ending L value
      time(k)=UT_end      !.. set the ending time
      NDint=1             !.. count # of values

      !.. loop to trace back from desired ending time to find the L values back
      !.. to the beginning of time
  10  CONTINUE
          !.. Interpolate the velocity
          Vint=FTAE(1,NWTIM,UT,Vion,time(k))   !... Get velocity at the closest time
          !..if(time(k) > UT(i)) i=i+1 

          !.. scale Vion to the equator
          !.. LSCAL is a scaling factor to reduce ExB drift with L_shell
          LSCAL=1.0
          IF(L(k).GT.1.05) LSCAL=EXP(LVSCAL*(L(k)-1.05))
          Veq(k) = Vel_Convert(L(k),Vint)*LSCAL
          Vsave(k)=Vint        !.. Save the input velocity for testing

          !.. Move the flux tube back to where it came from using Veq=DL/Dt
          L(k-1) = L(k)-Veq(k)*(time_step)*3600/6370000 

          !.. Make sure apex altitude > 190 
          IF(L(k-1).LT.(6370.+190.)/6370.) L(k-1)=(6370.+190.)/6370.

          IF(L(k-1).GT.PCO_MAX) PCO_MAX=L(k-1)     !.. find maximum L value
          IF(L(k-1).LT.PCO_MIN) PCO_MIN=L(k-1)     !.. find minimum L value

          NDint=NDint+1             !.. count # of values
          kmin=k-1                  !.. minimum value of k

          !.. Go to previous time
          k=k-1 
          time(k)=time(k+1)-time_step 

      IF(( k.GT.1) .AND. (time(k) > UT(1))) GOTO 10 !.. end of loop

      !... Move the data to start at first array position in MEDATA
      DO i=1,NDint
        j=kmin+i-1
        MEDATA(i,1)=time(j)
        MEDATA(i,2)=L(j)
        MEDATA(i,3)=Veq(j)
        MEDATA(i,4)=Vsave(j)
      ENDDO
        !.. Add extra time at the end to help reach end point
        NDint=NDint+1
        MEDATA(NDint,1)=time(NDint-1)+2
        MEDATA(NDint,2)=L(NDint-1)
        MEDATA(NDint,3)=Veq(NDint-1)
        MEDATA(NDint,4)=Vsave(NDint-1)
      RETURN
      END

C::::::::::::::::::::::::::: Vel_Convert:::::::::::::::::::::::::::::::::::
C..FUNCTION    : Vel_Convert 
C..DESCRIPTION : This function is used to scale Vion(ionospheric velocity) 
C.. to Veq(velocity at the equatorial plane) using a multiplication factor 
C.. from Mikhailov and Forster, 1998 JASTP paper. The formula is:
C.. Veq=Vion*((4.0-3.0/L)**0.5)*(L**1.5), using (cos(theta))**2 = 1/L.
C.. For programming simplicity, this formula is rewritten as: 
C.. Veq=Vion*L*((4.0*L-3.0)**0.5) 
C.. The formula is modified so that the scaling factor is unity in the ionosphere   
C.. AUTHOR      : Jinhua Liao in January, 1999, Modified for FLIP by 
C.. P. Richards Dcember 2008
C***********************************************************
      REAL FUNCTION Vel_Convert(L,    !.. L_shell for velocity calculation
     >                       Vion)    !.. ionospheric velocity 

      IMPLICIT NONE 
      REAL L,Vion
      REAL LNORM          !.. L normalization factor
      REAL ZNORM          !.. Normalization altitude
      DATA ZNORM/300.0/   !.. Normalization altitude

      !.. Factor to make Vel_convert scaling factor unity in the ionosphere at ZNORM
      LNORM=((6370.+ZNORM)/6370.)*sqrt(4.0*((6370.+ZNORM)/6370.)-3.0)

      !.. return the velocity at equatorial plane
      Vel_Convert=(Vion*L*sqrt(4*L-3))/LNORM
      RETURN
      END
